NASA TECHNICAL NOTE 


NASA TN D-7305 



AN EXTENSION OF THE QZ ALGORITHM 
FOR SOLVING THE GENERALIZED 
MATRIX EIGENVALUE PROBLEM 


by Robert C. Ward 

Langley Research Center 
Hampton, Va. 23665 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHIN6T0N, D. C. • JULT 1973 



1. Report No. 2. Government Accession No. 

NASA TN D-7305 

3. Recipient's Catalog No. 

4. Title and Subtitle 

AN EXTENSION OF THE QZ ALGORITHM FOR SOLVING THE 
GENERALIZED MATRIX EIGENVALUE PROBLEM 

5. Report Date 

July 1973 

6. Performing Organization Code 

7. Author(s) 

Robert G. Ward 

8. Performing Organization Report No. 

L-8899 

10. Work Unit No. 

501-06-01-11 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, Va, 23665 

11. Contract or Grant No, 

13, Type of Report and Period Covered 

Technical Note 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 


15. Supplementary Notes 

The basic information presented herein is a part of a thesis which will be offered in 
partial fulfillment of the requirements for the degree of Doctor of Philosophy in Applied 
Mathematics, University of Virginia, Charlottesville, Virginia. 

16. Abstract 


An algorithm called the combination shift QZ algorithm is presented for solving the 
generalized matrix eigenvalue problem. This new algorithm is an extension of Moler and 
Stewards QZ algorithm with some added features for saving time and operations. Also, 
some additional properties of the QR algorithm which were not practical to implement in 
the QZ algorithm can be generalized with the combination shift QZ algorithm. Numerous 
test cases are presented to give practical application tests for the algorithm. Based on 
the results presented in this paper, this algorithm should be preferred over existing 
algorithms which attempt to solve the class of generalized eigenproblems where both 
matrices are singular or nearly singular. 


Li 7. Key Words ^Suggested by Author{s)) 18. Distribution Statement 


Generalized eigenvalue problem 

Numerical analysis 

Eigenvalues 

Algorithms 

QZ algorithms 

Unclassified - 

- Unlimited 


19. Security Qassif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price* 

Unclassified 

Unclassified 

54 

$3.00 


For sale by the National Technical Information Service, Springfield, Virginia 22151 
























CONTENTS 


Page 

SUMMARY 1 

INTRODUCTION ........ ... 1 

SYMBOLS 3 

QZ ALGORITHM 4 

Major Steps of QZ Algorithm 5 

Step 2 of QZ Algorithm 6 

COMBINATION SHIFT QZ ALGORITHM 11 

Single Shift Implicit QZ Iteration 11 

Consecutive Small Subdiagonals . 15 

Step 2 of Combination Shift QZ Algorithm 19 

THEORETICAL COMPARISON OF THE COMBINATION SHIFT QZ 

AND THE QZ ITERATIONS 20 

NUMERICAL RESULTS 22 

CONCLUDING REMARKS 27 

APPENDIX - TEST CASES 28 

REFERENCES 44 

TABLES ............................... 45 


iii 



AN EXTENSION OF THE QZ ALGORITHM FOR SOLVING THE 
GENERALIZED MATRIX EIGENVALUE PROBLEM* 

By Robert C . Ward 
Langley Research Center 

SUMMARY 

An algorithm called the combination shift QZ algorithm is presented for solving the 
generalized matrix eigenvalue problem . This new algorithm is an extension of Moler and 
Stewart's QZ algorithm with some added features for saving time and operations. Also, 
some additional properties of the QR algorithm which were not practical to implement in 
the QZ algorithm can be generalized with the combination shift QZ algorithm. Numerous 
test cases are presented to give practical application tests for the algorithm . Based on 
the results presented in this paper, this algorithm should be preferred over existing 
algorithms which attempt to solve the class of generalized eigenproblems where both 
matrices are singular or nearly singular. 

INTRODUCTION 

There are numerous problems which occur frequently in the physical sciences that 
require solving the generalized eigenvalue problem 


Ax = XBx (1) 

for X and x where A and B are n Xn real matrices, X is a scalar, and x is 
a n X 1 vector. To mention one example, it is well known (see Lancaster (ref. 1)) that 
the equations of motion for many mechanical and electrical systems may be written in the 
matrix form 

Ap + Bp + Cp = f (2) 

where A, B, and C are n x n real matrices, and p and f are time -dependent 
n X 1 vectors . If a system with no damping (B = 0) , no forcing function (f = 0) , and 
solutions of the form p(t) = e^^ (sinusoidal solutions) where x is independent of 
time is considered, then in terms of X and x, equation (2) becomes 

*The basic information presented herein is a part of a thesis which will be offered 
in partial fulfillment of the requirements for the degree of Doctor of Philosophy in 
Applied Mathematics, University of Virginia, Charlottesville, Virginia. 
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Cx = -x2ax (3) 

which is the generalized eigenvalue problem in If the system is damped, the result- 

ing equation can be transformed by the equations 

Ax = y (4) 


AAy + By + Cx = 0 

into the block 2n x 2n matrix generalized eigenvalue problem 



(5) 


( 6 ) 


When solving the generalized eigenvalue problem, the roles of A and B can be 
reversed by solving for the reciprocals of the eigenvalues. That is, one could solve the 
problem 

Bx = ja Ax (7) 

where the eigenvalues A of the original problem (eq, (1)) are given by 

^ = ( 8 ) 

An infinite eigenvalue of equation (1) is defined as a corresponding zero eigenvalue of 
equation (7) . 

In many application cases, A and B of equation (1) have some special properties 
which determine the type of eigenvalues present and influence the selection of an algorithm 
for solving the problem. First, consider the cases when B is well- conditioned with 
respect to inversion. When A is symmetric and B is positive definite as quite often 
occurs in physical applications, all the eigenvalues are real; Martin and Wilkinson (ref . 2) 
describe an algorithm which reduces this problem to the standard symmetric eigenvalue 
problem Pz = Az. In addition, if A and B have band structure, Crawford (ref. 3) 
describes a modification to Martin and Wilkinson's algorithm to take advantage of the 
band matrices. Also, an algorithm by Golub, Underwood, and Wilkinson (ref . 4) using 
Lanczos method solves the band problem. If A and B do not have special proper- 
ties, the generalized eigenvalue problem can be solved by solving the standard problem 
B"^Ax = Ax. Normally, one would not want to form B'^A when A and B have spe- 
cial properties usable by algorithms. For example, in the preceding case where A was 
symmetric and B was positive definite, the problem was transformed into a standard 



symmetric eigenvalue problem which is faster and numerically more stable than the non- 
symmetric problem B" ^Ax = Ax. 

Now, consider the more complicated cases when B is ill-conditioned. When A 
is symmetric and B is positive semidefinite or positive definite but ill-conditioned, 

Fix and Heiberger (ref . 5) describe an algorithm for solving this problem which depends 
on determining the rank of several submatrices. This case is interesting in that the 
spectrum consists of both stable and unstable real eigenvalues and there exists the possi- 
bility of every scalar being an eigenvalue . Unstable eigenvalues are those which are sen- 
sitive to small changes in the matrix elements of A and B and thus cannot be computed 
accurately by a normal computational procedure. In general, rank determination is a dif- 
ficult problem to solve on a computer; thus. Fix and Heiberger* s algorithm runs into diffi- 
culty when there is not a clear separation between the stable and unstable eigenvalues. 

For a general ill-conditioned B matrix, Peters and Wilkinson (ref, 6) describe an algo- 
rithm which again depends on rank determination. Moler and Stewart (ref. 7) describe an 
algorithm which solves the generalized eigenvalue problem for arbitrary real matrices A 
and B by use of unitary transformations. This generalization of the double shift QR 
(ref. 8) is called the QZ algorithm and is particularly effective for the cases when B is 
singular or nearly singular . The algorithm presented in this paper is a combination of 
the QZ and a generalization of the single shift implicit QR and is referred to as the com- 
bination shift QZ algorithm. This algorithm is also effective on the singular or nearly 
singular B cases and is particularly effective on cases where a large number of real 
eigenvalues are expected, such as A symmetric and B positive semidefinite. 

SYMBOLS 


A,B,C 

a’,b’,c' 


n by n matrices 

n by n matrices, next iterate of A, B, and C, respectively 


A,B 


lower right 2 by 2 submatrices of A and B, respectively 


a|j,b|j i,j elements in matrix A and B, respectively 


3-10»^20>^30»l 
S2,u,r 


scalars 


G n by n matrix related to C 

Dj^,D| block diagonal matrices 
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p,x 


I 

k,Z,m,q 

n 


Q,Z 


Q’,Z’,Z" 

Qi,Zi 


H>yi 




^0 

V2’Vi,V2 


n by 1 vectors 
identity matrix 

index arguments for derivation of test eases 
size of matrices A and B 
n by n orthogonal matrices 
n by n orthogonal matrices 
ith n by n orthogonal matrix in a sequence 
n by 1 vectors associated with eigenvalue Xj 
scalars 

basic machine roundoff error 
scalars, usually small numbers 


cr shift in single shift implicit QZ iteration 

shifts in double shift QZ iteration 



the norm of [ ] 

Hermitian of [ j 
transpose of [ j 

QZ ALGORITHM 


Since a detailed description of the QZ algorithm is given in reference 7, only a brief 
summary will be given here for completeness. The algorithm is an iterative method for 
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computing the decomposition guaranteed in the following theorem from reference 9: 

Theorem : There are unitary matrices Q and Z so that QAZ and QBZ are both upper 
triangular. 

If the decomposition can be accomplished, then the eigenvalues and eigenvectors 
are easily extracted by 


Xi = Zy. (10) 

where and jSj are the diagonal elements of the QAZ and QBZ matrices, respec- 
tively, and y. are the eigenvectors of the triangular system QAZy^ = XjQBZy^. 

Major Steps of QZ Algorithm 

There are four major steps in the algorithm. 

(1) Reduce A to upper Hessenberg form and at the same time reduce B to upper 
triangular form. 

(2) Use a generalization of the double shift QR to put A in quasi-triangular form 
(upper block triangular form with 1 x 1 or 2 x 2 diagonal blocks) while keeping B in 
upper triangular form. 

(3) Reduce A to upper triangular form and keep B in upper triangular form . 

(4) Find the eigenvectors of the triangular system and back transform them to the 
original problem . 

Step 1 is initiated by transforming B into upper triangular form by premultiplying 
by a unitary matrix, denoted Q*, made up of a sequence of Householder reflections. Then 
Q’A is put in upper Hessenberg form by annihilating one element at a time in the order 
given below on a 5 x 5 example: 
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After each annihilation, for example by Q^, the current B matrix is put back in triangu- 
lar form by a postmultiplication with a Householder reflection, for example Zj. This 
reflection does not affect the zeros introduced in the current A matrix. Thus, in the 
5x5 matrix example, QAZ and QBZ would be upper Hessenberg and upper triangular, 
respectively, where Q = QgQ 5 Q 4 Q 3 Q 2 QiQ' and Z = Z 2 Z 2 Z 3 Z 4 Z 5 Zg. 

Step 2 is the generalization of the double shift QR algorithm applied to the standard 
eigenvalue problem AB"^x = Ax without forming B“^> Since this step is the iterative 
step that the combination shift QZ algorithm alters, a more detailed description of this 
iteration is given later. For a complete description, the interested reader should con- 
sult reference 7. 

Step 3 involves one QZ transformation, that is, one premultiplication by a unitary 
matrix Q and one postmultiplication by a unitary matrix Z on the quasi-triangular 
A and triangular B to put both matrices in triangular form. If a 2 x 2 diagonal block 
of A corresponds to a complex conjugate eigenvalue pair, complex arithmetic will be 
required in this step. 

Step 4 is accomplished by solving the reduced triangular problem for its eigenvec- 
tors by a back- substitution process similar to the one used by Peters and Wilkinson 
(ref. 10) in the procedure "hqr 2.” The Z transformations (postmultiplication matri- 
ces) are accumulated and applied to the eigenvectors of the reduced system to obtain the 
eigenvectors of the original system. Recently, Kaufman (ref. 11) has pointed out that it 
is advantageous to solve the transposed problem for the left eigenvectors and thus accu- 
mulate the Q transformations for back transforming the triangular system vectors. 

This advantage can be seen in the discussion of the second step in the QZ algorithm . 

Step 2 of QZ Algorithm 

The iteration is motivated by assuming that B is nonsingular and by examining 
the double shift QR algorithm for C = AB"1 . Recall that A is upper Hessenberg and 
B is upper triangular as a result of step 1; thus, G is also upper Hessenberg. 

Suppose one iteration of the double shift QR with shifts and Og is applied to 
C . Then a unitary matrix Q is found that makes the matrix QCQ® upper Hessenberg, 
where Q^i denotes the Hermitian of Q, and the first row of Q is the first row of the 
orthogonal matrix which annihilates all the elements but the first in the first column of 
- Oglj. The next iterate C' is then defined as 

C’ = QCQH , (11) 

Consider what happens if a special form of the identity, that is, ZZ^ where Z 
is unitary is inserted in the matrix equation (1) to give 
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Then 


Az(z*^) = XBz(zHx) 


( 12 ) 


C = AZ(BZ)-1 = AZZ%"^ = AB-1 = G (13) 

Since the unitary matrix Z does not change C or C' , Z could be used effectively to 
keep A upper Hessenberg and B upper triangular during the iteration and not destroy 
the zeros introduced in step 1, 

Suppose the matrix Q in equation (11) is known and the matrix Z which keeps 
the current A and B in the proper form is known; Define A' and B' by 

A' = QAZ (14) 

B’=QBZ (15) 

and then 


A’B'-l = (QAZ)(QBZ)'l = QAZZ^B-Iq^ = QAB'^QH = C' (16) 


Thus, if Q and Z can be determined without forming B"^, then the next iterates, A' 
and B', can be determined by equations (14) and (15). 


The QZ iteration must then do two things. First, determine the correct first row 
of Q. Second, determine Q and Z so that Q retains the correct first row, QAZ is 
upper Hessenberg, and QBZ is upper triangular. 


As mentioned earlier, the first row of Q is the first row of the orthogonal matrix 
which annihilates all the elements but the first in the first column of ^AB"1 - a^lj^AB"! - 
a 2 l). Since A is upper Hessenberg and B is upper triangular, the first column of 
^AB"1 - a^lj^AB"^ - o- 2 l^ is completely determined by Oj, CT 2 , and the first two columns 
of AB"1. Only the nonsingularity of the upper 2 by 2 submatrix of B, that is, bj^j and 
b 22 nonzero, is required to find these columns. This is really no restriction since Moler 
and Stewart show how to handle the case of the singular or nearly singular submatrix. 


The first column of ^AB~1 - CTj^I^^AB"^ - is called "the fictitious zeroth column" 
of AB”1 and is easily computed. The first column of AB"1 has two nonzero elements 
and the second column has three, The equations for these elements are 


(ab-i)u = \ 


11 

11 


(ab-i),j = 


21 “ b 


£l 

11 


(17) 

(18) 
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( 19 ) 


(ab-1) , ^12 ^11^12 

t»22 ^ 1^22 


rAR-1) -^22 ^2l'"l2 

^22 bgg ^11^22 


( 20 ) 




22 


( 21 ) 


Since the rate of convergence for the QR algorithm is determined by the ratio of consecu 
tive distinct eigenvalues, shifts are employed to make this ratio as small as possible. 
The shifts and Og are chosen to be the two zeros of the equation 

det(A - ob) = 0 

where 


( 22 ) 


A = 


a 


tt-l,n-l 

a - 
n-l,n 

a < 

n,n- 1 

^n,n 

*n-l,n-l 

Vl,n 

0 

‘’n.n 


B = 


and n is the order of the current A and B matrices. These shifts are not explicitly 
computed, but techniques similar to those used in "hqr2" (ref. 10) are applied instead. 

By denoting n - 1 by m, the three nonzero elements of the zeroth column of AB"1 
are computed by the following formulas: 


aio = 


*mm 

^mm 


in 

^22 



^llA’^nn 



'■21y 


(23) 


a 


20 


_ f^22 ^11^ 


^22 


^30 =b 


_ ^32 


22 


^^2lV^12^ 

^hl\^22j 


mm 

^mm 



(24) 

(25) 
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Note that Ej^q, a2o» and a^Q are all real even if the roots (shifts) of equation ( 22 ) are 
complex. This is an important feature of the QZ algorithm. 

After a|Q, a2Q, and a0Q have been computed, the iteration involves only premul- 
tiplying and postmultiplying the A and B matrices by Householder transformations, 
the proper elements being annihilated each time. Since there are only three nonzero ele- 
ments in the "zeroth" column, only the first three elements of the first row of Q are 
nonzero. Being illustrated on 6 by 6 matrices, A and B have the following form after 
applying the Householder transformation for annihilating a2Q and agQ: 


QjA QjB 



The algorithm must now reduce Qj^A to upper Hessenberg and Qj^B to upper triangu- 
lar without affecting the first row. This is done by postmultiplying by the Householder 
transformation which annihilates the elements denoted by superscript 1 and then by 
the Householder transformation z'][ which annihilates the element in position denoted by 
superscript 2 , Letting Zj^ - Z^Z'^ yields the following forms: 


QjAZ^ Q^BZ^ 
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Now, annihilating the elements denoted by superscript 1 with a Householder transforma- 
tion Q 2 yields 


QgQjAZ^ QgQiBZj 



Postmultiplying by Z 2 » a product of Householder transformations annihilating the ele- 
ments in the given order, reduces the current B matrix to triangular form. Then pre- 
multiplying by Qg annihilates the nonzero elements outside the Hessenberg form in the 
second column of Q 2 Qj^AZ^Z 2 * This procedure continues until all the unwanted nonzero 
elements are pushed down to the lower right-hand corner and a Hessenberg matrix A' 
and an upper triangular matrix B* remain. 

By letting the elements of the current transformed A and B matrices be denoted 
by ay and by, respectively, the iteration can be summarized by the following outline: 

(1) Compute aj^Q, agg, and agg by equations (23), (24), and (25) 

(2) For k = 1, 2, . . ., n - 2, 

(a) Determine to annihilate 3.^^^ i and k- 1 

(b) Determine Z^ to annihilate hjj ^.2 and k 

(c) Determine z'^ to annihilate bj^^j 

(3) Determine to annihilate a^^ j ^_2 

(4) Determine Z„ ■> to annihilate b„ „ ■< 

' ' n-i n,n-l 

The operation count (only operations of the highest order of n are counted) for one 
double shift iteration is 13n2 multiplications, 13n2 additions, and 3n square roots. 

If the eigenvectors are required, then the Z matrices must be accumulated, which adds 
8n2 more multiplications and 8n2 more additions per iteration. K the transposed prob- 
lem is solved, then the Q matrices are acciimulated and add 5n2 more multiplications 
arid 5n2 more additions per iteration instead of 8n2. 
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COMBINATION SHIFT QZ ALGORITHM 

The combination shift QZ algorithm is basically the QZ algorithm with two improve- 
ments which take advantage of some opportunities for saving time and operations. Steps 1, 
3, and 4 are not altered by the new algorithm. Thus, the iterative step which is the heart 
of the algorithm is the only step affected. 

Step 2 was the generalization of the double shift QR. The double shift QR is used 
to solve the standard eigenvalue problem AB'^x = Ax because of the possibility of com- 
plex conjugate shifts and Og. If the shifts are complex, the double shift version 
allows the continuation of the use of real arithmetic,, as pointed out earlier. However, 
if the shifts are real, this feature of the double shift version is no longer an advantage. 
Thus, a generalization of the single shift implicit QR algorithm might have some advan- 
tages when real shifts are encountered. After a discussion of this generalization and 
one of its properties that can be utilized, the second step of the new algorithm will be 
explained. 


Single Shift Implicit QZ Iteration 

Similar to the double shift generalization, the single shift implicit QZ iteration is 
motivated by assuming B is nonsingular and by exaniining the single shift implicit QR 
algorithm for C = AB"^, Recall that C is an upper Hessenberg matrix because of 
step 1. 

Suppose one iteration of the single shift implicit QR algorithm with shift a is 
applied to C . Then a matrix Q is found so that the matrix QCQ® is upper 
Hessenberg and the first row of Q is the first row of the orthogonal matrix which 
annihilates all the elements except the first in the first column of the matrix (C - al). 

The next iterate C’ is then defined as 

C'=QCQH (26) 

As was the case in the double shift, A and B can be postmultiplied by a unitary 
matrix Z without altering C or C' . Thus again this technique can be used to keep 
A upper Hessenberg and B upper triangular during the iteration. Also by finding the 
proper matrices Q and Z, the next iterates A' and B' can be found without explic- 
itly forming B"^ by using the equations 

A' = QAZ (27) 

B’ = QBZ (28) 
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Thus, to perform one single shift implicit QZ iteration, the algorithm must do two 
things. First, determine the correct first row of Q. Second, determine Q and Z so 
that Q retains the correct first row, QAZ is upper Hessenberg, and QBZ is upper 
triangular. 

The first row of Q is the first row of the orthogonal matrix which annihilates all 
the elements except the first in the first column of (aB“ 1 - al). Since A is upper 
Hessenberg and B is upper triangular, the first column of (AB"1 - ol) is completely 
determined by o, a 2 i, and bj^^. In fact, the first column of (AB“ ^ - alj, called 

the fictitious zeroth column, has as its first two elements 


11 


(29) 


( 30 ) 

with the remaining elements all zero. Thus, a nonzero b^j^ is the only requirement on 
the nonsingularity of B. If bj^^ is zero, then a deflation can be carried out to reduce 
the order of the working matrices A and B, This procedure will be discussed later. 

The second part is very similar to the second part of the double shift iteration. 
Premultiplication and postmultiplication by Householder transformations are alterna- 
tively used to annihilate the proper elements to reduce A to upper Hessenberg form 
and B to upper triangular form without affecting the first row. Let Qj be the uni- 
tary matrix which annihilates a 2 Q| that is, the first row of Q| is the desired first 
row for Q. Then Qj^A and Qj^B have the following form on a 5 by 5 example: 


Q^A 


QlB 

X X X X X 


X X X X X 

X X X X X 


X X X X X 

0 X X X X 


0 0 X X X 

0 0 X X X 


0 0 0 X X 

0 0 0 X X 


0 0 0 0 X 


Now, Q^B must be returned to upper triangular form while Qj^A is kept in upper 
Hessenberg form. Postmultiplying by to annihilate the element in the b 2 j^ posi- 
tion yields 
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Premultiplying by Qg to annihilate the element in the sl^i position yields 


QgQjAZj QgQiBZj 



Postmultiplying by Zg to annihilate bgg yields 


Q„Q,AZ^Z„ Q„Q,BZ,Z„ 

2112 ^21 1 2 



The process continues with Qg annihilating a^2» Zg annihilating b^g, annihi- 
lating agg, and finishing with Z^ annihilating bg^. This procedure yields a unitary 
matrix Q = Q4Q3Q2^1’ ^ unitary matrix Z = Z^ZgZgZ^, an upper Hessenberg matrix 
A* = QAZ, and an upper triangular matrix B’ = QBZ so that the first row of Q is the 
first row of Qj as required. 
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By letting the elements of the current transformed A and B matrices be denoted 
by ay and by , respectively, the single shift implicit QZ iteration can be summarized 
by the following outline: 

(1) Compute aj^Q and a 2 Q by equations (29) and (30) 

(2) For k = 1, 2, . , n - 1, 

(a) Determine Qj^ to annihilate l 

(b) Determine Zjj. to annihilate b^^^ 

The operation count for one iteration is 6n2 multiplications, 6n2 additions, and 
2n square roots. K the eigenvectors are required and the Z matrices are accumulated, 
this procedure adds 3n2 more multiplications and 3n^ more additions per iteration. 
There is no advantage in solving the transposed problem and accumulating the Q matri- 
ces for this iteration since it would require the same number of operations. This condi- 
tion is due to only one Z transformation being required to annihilate the B matrix 
elements as they become nonzero in this iteration whereas two Z transformations are 
required in the double shift QZ iteration. 

As previously noted, the iteration runs into trouble if bi ^ is zero or negligible. 

^A matrix element is defined as negligible if the element in absolute value is less than 
6 q times the norm of the matrix. The infinity matrix norm is used in both the QZ and 
the combination shift QZ algorithms, The error €q Is the basic machine roundoff 
error. (See Wilkinson (ref . 12) .)) The solution to this problem is to deflate from the 
top as is done in the double shift QZ . If bj^j is negligible, it may be set equal to zero 
without difficulty since unitary transformations are being used. For a 4 by 4 example, 

A and B then have the form 


A 


B 

X X X X 


0 X X X 

X X X X 


0 XXX 

0 X X X 


0 0 X X 

0 0 X X 


0 0 0 X 


Thus a unitary Q can be used to annihilate agj^ without affecting the form of B. This 
procedure gives a zero subdiagonal element of A and reduces the problem. If bjj is 
not quite small enough to be set equal to zero, it will cause the shift to be felt only weakly 
since the division by bj^j will overshadow o in the equation for ajo and a20- How- 
ever, the iteration can still be used profitably in converging to a large eigenvalue at the 
top. 
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Both the double shift QZ and the single shift QZ have the property of reducing the 
problem when a negligible subdiagonal element of A is encoimtered. On a 6 by 6 exam- 
ple, if a .^2 negligible, it can be set equal to zero and the eigenvalues of the full matrix 
problem can be found by solving for the eigenvalues of the following two diagonal systems: 


A B 



But because of the simplicity of the single shift iteration, one more generalization of the 
basic QR algorithm can be applied. The property of reducing the number of transfor- 
mations by detecting two consecuttve small subdiagonal elements of A can now be 
generalized. 

Consecutive Small Subdiagonals 

Suppose Uj. j._j| and a^,^! j. are both "small," but not negligible. One would like 
to develop a test similar to the QR algorithm that would allow the iteration to start at col- 
umn r, that is, the Qjj matrices would affect only rows r and below. Let the form of 
A and B on a 6 by 6 example for r = 3 be 


A B 



where ej^(a32) and 62(^43) "small” in some sense. 
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Suppose the iteration starts at column 3 instead of column 1. Then by making the 
assumption that b 33 is nonzero or nonnegligible, one has 


^33 

^33 


a 


(31) 


^20 - bgg 


(32) 


A discussion will be given later on what can be done if this assumption is false. Let the 
transformation Qg which annihilates agQ be denoted as 

Qg = I - p(o!p’^) (33) 

where the vector p and scalar a are found by the following set of equations: 

pT = (0, 0, l,u, 0, 0) (34) 


S2 = 



H- 



(35) 


^ ^20 
Ufo ± S 


(36) 


a 


2 

1 + u2 


(37) 


The sign of S in equation (36) is chosen so that S and aj^Q have the same sign. 
Thus, QgA and QgB would have the forms 

Q3A QgB 


X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


0 

X 

X 

X 

X 

X 

0 

^1 

X 

X 

X 

X 


0 

0 

X 

X 

X 

X 

0 

^2 

X 

X 

X 

X 


0 

0 

X 

X 

X 

X 

0 

0 

0 

X 

X 

X 


0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

X 


0 

0 

0 

0 

0 

X 
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where 


Vi = - 


2^1 
1 + u2 


(38) 


^2 = - 


2ej^u 
1 + u2 


(39) 


H 172 i® negligible, then one is justified in starting the iteration with the third col- 
umn. One would set ryg equal to zero, apply Zg to annihilate b^g and introduce a 
nonzero in a^g, apply to annihilate a^g and introduce a nonzero in bg^, and so 
forth. Hence, an easily computable bound on jTyg j must be obtained. From equa- 
tions (36), (35), and the sign selection, one finds 


1 

ju ] = 

^20 

< 

^20 


^20 

^10 ^ ® 


^10 ± \/^ 


2 a 10 


(40) 


By using equations (39) and (40), 


^2 


2ej^u 

< 

Op n 

< 

^1^20 


l+u2 




^10 


(41) 


From equations (31), (32), and (41), 


€ i 62 


^32^43 

^33 “ ^^33 


^33 '■ °"^33 


(42) 


Equations (40), (41), and (42) require a^Q to be nonzero. It is conceivable for 
aj^Q to be zero. If this occurs, rj 2 would be equal to ±€j^ and would be nonnegligible. 
Therefore, the following test inserting the general index r, is used for determining the 
negligibility of tj 2 instead of equation (42) : 


a. 


r,r-l^r+l,r ~ 


«-rr 


ab 


rr 


€o||A| 


(43) 


If this inequality holds, the iteration can start in column r instead of column 1. In this 
specific example, the iteration can start in column 3 instead of column 1. 

Now, consider equation (38) and try to derive an expression for when Tjg is 
negligible 
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(44) 


77j = ej- 


2ei - 6 j^u 2 + 2€ju2 


1 + u2 1 + u2 


Thus, one has 


n. = - 


6i + 


2ej^u 


1 + 


u = -e, - 


1- V 


From equation (35), 


= |^20j 

which yields from equation (36) and the sign selection 


(45) 


(46) 


^20 

< 

*20 

aio ± S 


|Ho| + |^2of 


(47) 


Since is negligible, one now has an expression for r}^ 

7?i = -6^ = -agg (48) 

which just involves changing signs. 

There is only one difficulty to clear up. The fact that bgg is not negligible was 
used even though the final result that is tested (eq. (43)) does not require this. K bgg 
is negligible, then one may set bgg equal to zero and try to perform a reduction. After 
setting bgg equal to zero, A and B have the forms 


A B 


X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


0 

X 

X 

X 

X 

X 

0 


X 

X 

X 

X 


0 

0 

0 

X 

X 

X 

0 

0 

^2 

X 

X 

X 


0 

0 

0 

X 

X 

X 

0 

0 

0 

X 

X 

X 


0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

X 


0 

0 

0 

0 

0 

X 
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It is desirable to annihilate €2 by a Householder transformation to reduce the problem. 
By setting 


10 " ^33 

(49) 

20 = ^2 

(50) 


and going through the analysis of equations (33) to (41), in this case obeys the fol- 
lowing inequality: 



^ 1^2 

^33 


(51) 


If ^2 negligible, 77 ^ can be found through the same analysis as equations (44) to (48) 
with the same result 



(52) 


Thus, if 


^ 1^2 


'33 


is negligible, the sign on 


problem reduced. 


can be changed, €2 annihilated, and the 


In summary, the algorithm tests the subdiagonal elements a^. and a^,^j ^ 

(r = n - l,n-2, . . ., 2) each iteration for the validity of equation (43) if bj.^. is nonneg- 
ligible, and equation (43) with replaced by japj.| if bj.j. is negligible. 

Suppose equation (43) is valid for subdiagonal elements a^^ and a^^j^ The algo- 
rithm then proceeds according to the negligibility of b^^. If bji is not negligible, the 
iteration starts a column i. If bj| is negligible, the matrix problem is reduced into 
two smaller matrix problems . 


Step 2 of Combination Shift QZ Algorithm 

The second step of the combination shift QZ algorithm can now be stated; that is, 
reduce A to quasi- triangular form while keeping B in upper 
triangular form by using a combination of the double shift QZ 
and the single shift implicit QZ. 

The type of iteration used to determine the next iterates depends on the type of 
shifts computed. By using equation (22), the algorithm determines whether the shifts 
are real or complex. If they are complex, a double shift QZ iteration is performed as 
explained earlier. If they are real, a single shift implicit QZ iteration is performed by 
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using the shift which is closest to the value where and are the n,n 

°nn 

elements of the current transformed A and B matrices, respectively. 

The reasoning behind the selection of this step is obvious. K the shifts are complex, 
then the calculation remains in real arithmetic by the double shift and is probably con- 
verging to complex eigenvalues, If the shifts are real, then the Iteration is probably con- 

^nii 

verging to at least one real eigenvalue which will emerge as t^. 

“nn 

Since the double shift QZ may be used, it is still an advantage to solve the trans- 
posed problem and accumulate the Q transformations. This method will insure only 
having to accumulate one Householder reflection per iteration. 

THEORETICAL COMPARISON OF THE COMBINATION SHIFT QZ 

AND THE QZ ITERATIONS 

One important and interesting form of comparison is operation count (only multi- 
plications and divisions of the highest order of n are stated here). The operation count 
for one double shift QZ iteration is 13n2, whereas the operation count for one combina- 
tion shift QZ iteration depends on the type of shifts encountered. If the shifts are com- 
plex, then the count is the same as that of the double shift. If the shifts are real, then 
one shift is performed and the iteration requires 6n2 operations. (One should note that 
there are a few more logical statements and multiplications in the combination shift itera- 
tion because of the shift-type determination, but these are of order unity per iteration.) 

To emphasize what happens in the real case, suppose the two algorithms are con- 
verging to a real eigenvalue and the shifts in both the double shift QZ and the combination 
shift QZ are real. The combination shift QZ would iterate with 6n2 operations, obtain 
a new shift estimate, iterate again with 6n2 operations, obtain a new shift estimate, and 
so on until convergence. The double shift QZ would use both shifts and iterate with 13n2 
operations, obtain two new shift estimates. Iterate again with 13n2 operations, and so 
on until convergence. Thus, the combination shift QZ can perform two iterations with a 
better shift estimate for the second iteration with n2 fewer operations than the double 
shift QZ. Consequently, when trying to converge to a real eigenvalue, it appears that the 
combination shift would save on the number of operations per two shifts as well as possi- 
bly on the number of iterations because of the improved shift estimates. When trying to 
converge to a complex eigenvalue pair, the shifts should be complex and the two algorithms 
should be roughly equivalent. One would then expect that the combination shift QZ itera- 
tion would be faster than the double shift QZ iteration when real eigenvalues are present, 
but would at worst (all complex eigenvalues) be only slightly slower. One would also 
expect that the savings due to the combination shift would be somewhat proportional to 
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the number of real eigenvalues found. These expectations are confirmed by the numeri- 
cal results presented later . 

The numerical tests also indicate that on large matrices, for example of order 50, 
time and iterations are generally gained or lost in computing the first 50 percent of the 
eigenvalues. There are two reasons for this gain or loss. First, the size of the matrix 
is reduced as an eigenvalue or an eigenvalue pair is found. Thus, there are more opera- 
tions associated with finding the earlier eigenvalues, and an iteration saved at the begin- 
ning is worth more than an iteration saved at the end. Second, the earlier iterations help 
the later iterations by orienting the eigenvalues in approximate order and giving better 
estimates for the shifts. (This is also the reason that extra early iterations are not com- 
pletely wasted.) In large matrices, the average rate of determining eigenvalues is gen- 
erally one eigenvalue per two shift iterations (one double shift or two single shifts) or 
better after the first 50 percent of the eigenvalues are found. It would then appear that 
the combination shift QZ algorithm would benefit from all the real eigenvalues being in 
a position to be found first. This advantage is somewhat offset by the property discussed 
in the next paragraph. 

Orientation of the eigenvalues plays a further role in determining the relative merits 
of the combination shift iteration. If a real eigenvalue would normally be found between 
two complex pairs of eigenvalues, the combination shift would operate more efficiently 
than the QZ since it has the capability of finding just one real eigenvalue. The double 
shift algorithm would have to disorder the eigenvalues in order to make use of both shifts 
in the iteration, or would have to perform an iteration with one shift which does not give 
immediate help in extracting the eigenvalues. To illustrate the point, consider the fol- 
lowing example: 

A B 


0 -3 -3 


10 -3 

1 1 -2 


Oil 

0 1 -2 


0 0 1 


The eigenvalues of Ax = XBx are -3, i ± i. The two shifts and used in the 
double shift are -2 and 0. The shift -2 can be used to good advantage in finding the eigen- 
value -3, but the shift 0 does not provide much help in converging to the eigenvalues. The 
combination shift would iterate by using a shift of -2 and would obtain a new shift for the 
next iteration which should be an even closer approximation to -3. This property of the 
- algorithms indicates that the combination shift should operate more efficiently than the 
QZ when the real eigenvalues are not bunched together. But when the real eigenvalues 
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are bunched, both shifts in the double shift could be used advantageously and the savings 
of the combination shift QZ is not as large. 

As previously mentioned, the simplicity of the single shift iteration allows one to 
check for consecutive small subdiagonals. Of course, the time and operations saved by 
utilizing this property is a function of the number of consecutive small subdiagonals 
detected by the algorithm and the position of these small elements along the subdiago- 
nal. In the 180 test cases presented in this paper, consecutive small subdiagonals were 
detected on the average of once every 8.7 single shift iterations. Since the double shift 
iteration does not have this capability without consuming considerably more time, this is 
a positive feature for the single shift iteration. 

The single shift iteration had to deal with the problem of a possible negligible bj j 
element. This problem was solved by a deflation from the top; that is, finding an infinite 
eigenvalue and reducing the order of the matrix problem by one, K bjj^ was small but 
not negligible, then an iteration was performed with a shift essentially equal to zero. 

This shift was useful in finding the large eigenvalue and thus deflating from the top, but 
did not provide much help In converging to the stable eigenvalue and deflating at the bot- 
tom. The double shift iteration has the same problem with bj| and the additional prob- 
lem of a negligible or almost negligible b 22 - The double shift iteration cannot take 
advantage of a negligible b 22 a-nd is forced to perforni an iteration with a shift essen- 
tially equal to zero when b 22 is negligible as well as almost negligible. Hence, the 
combination shift QZ algorithm would not be as likely to perform an iteration which does 
not help in converging to a stable eigenvalue as the QZ algorithm . 

Since both iterations involve only unitary transformations, they are both stable and 
well defined and, as expected, the accuracies of the two different iteration strategies are 
roughly equivalent. 


NUMERICAL RESULTS 


In order to determine the relative merits of the combination shift QZ algorithm, 
numerous test cases were run on the Control Data Corporation (CDC) 6600 computers 
at Langley Research Center . The results of these test cases were compared with the 
results of the same test cases by using Moler and Stewart's QZ algorithm. The QZ 
algorithm used was a FORTRAN computer code supplied by Dr. Cleve Moler of Univer- 
sity of New Mexico. For the combination shift QZ algorithm, only the subroutine of this 
code involving the iteration (that is, step 2 of the algorithm) was modified. The test 
cases were divided into six categories depending on the percentage of real eigenvalues 
possessed by the test case. The appendix gives the details on the generation of the 
matrices for these tests. Iteration times for the different test cases are given in 
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tables I to VI. Operation count and total algorithm time for category I are given in 
tables VII and VIII, respectively. 

The first category is one which is often found in the physical sciences; that is, the 
generalized eigenproblem with all real eigenvalues. For all the test cases in this cate- 
gory, A is a symmetric matrix and B is a symmetric, positive definite, nearly sin- 
gular matrix. For all the test cases in this category except 1-6, A is also nearly 
singular. Cases I-l, 1-2, 1-3, and 1-6 test the algorithms on problems which consist 
mainly of stable eigenvalues. The problem has two unstable eigenvalues, in test cases 
I-l, 1-2, and 1-3 and three unstable eigenvalues in test cases 1-6. Cases 1-4 and 1-5 test 
the algorithms on problems consisting of an approximiateiy equal number of stable and 
unstable eigenvalues. Most of the stable eigenvalues in these two cases are nearly zero. 
Table I reports the percentage of the QZ iteration time required by the combination shift 
QZ iteration. For the 30 test cases tried, the average time saved was over 35 percent 
and the deviation from this time was small. However, there was a definite trend toward 
larger savings on smaller matrices. To give an operation count comparison, a counter 
was inserted into both the combination shift QZ iteration and the QZ iteration to count all 
multiplications and divisions of order unity and above per iteration. Table VII reports 
the result of this comparison. As one would expect, results similar to table I are obtained 
with slightly larger percentages saved. Table Vni reports on the time comparisons of 
the complete algorithms. Since the algorithms are identical except for the iteration step, 
the results are also similar to those reported for the iteration time except with the per- 
centages closer to 100, as expected. 

The second category is the one with the test cases which have 80 percent real eigen- 
values. Table II reports the iteration time comparison for this category. The test 
cases n-1 have complex eigenvalues with larger magnitudes than the real eigenvalues. 

The test cases II-2 are just the opposite with the complex eigenvalues having smaller 
magnitudes than the real eigenvalues. Test cases II- 3 and II- 4 are problems which have 
each complex conjugate pair of eigenvalues isolated in magnitude; that is, if the eigenval- 
ues were ordered by magnitude, there would be at least one real eigenvalue between every 
complex conjugate pair. Test cases II- 5 and n-6 are problems which have the complex 
eigenvalues grouped together by magnitude; that is, if the eigenvalues were ordered by 
magnitude, there would be more than one complex conjugate pair between the real eigen- 
values. For the 30 test cases tried in this category, the average time saved in the itera- 
tion section of the algorithm is over 25 percent. Again, there is a trend for saving more 
on the smaller matrices. Also, the deviations are larger on the smaller matrices because 
of the larger effect of just one more or less iteration. 

The third category contains the test cases with 60 percent real eigenvalues . 

Table HI reports on the iteration time comparison between the QZ and the combination 
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shift QZ. Again, test cases III-l have complex eigenvalues with larger magnitudes, in-2 
has complex eigenvalues with smaller magnitudes, HI- 3 and HI- 4 have complex eigenval- 
ues isolated, and ni-5 and in-6 have complex eigenvalues grouped. For the 30 test cases 
tried in this category, the average iteration time saved is almost 30 percent. 

Table IV and table V report the results of the fourth and fifth category, respectively. 
The fourth category contains matrices with 40 percent of the eigenvalues real and the fifth 
category with 20 percent of the eigenvalues real. Again, test cases IV-1, V-1, IV-2, and 
V-2 have the same property as the respective cases in the two previous categories. Test 
cases IV-3, V-3, IV-4, V-4, IV-5, V-5, IV-6, and V-6 have their real eigenvalues ordered 
like the complex eigenvalues are ordered in the respective cases of categories II and III. 

In the fourth category, a savings of almost 20 percent is realized in the iteration time. 

The average iteration time saved in the fifth category is 13 percent. 

The sixth and last category is that of the all complex eigenvalue cases. Table VI 
reports the results of this category. Cases VI- 1 and VI-2 test the algorithms on prob- 
lems which consist mainly of stable purely imaginary eigenvalues. Cases VI- 3 and VI-4 
test the algorithms on problems which have eigenvalues with small and large imaginary 
parts relative to the real part and eigenvalues with real and imaginary parts of the same 
order of magnitude. Cases VI- 5 and VI- 6 test the algorithms on problems consisting 
of an approximately equal number of stable and unstable purely imaginary eigenvalues. 
Most of the stable eigenvalues in these two cases are nearly zero. An explanation for 
the savings in this category is that Some real shifts were encountered during the itera- 
tion and the combination shift QZ algorithm returned to the shift determination strategy 
quicker, and thus avoided as many real shifts as the QZ algorithm. Also, small real 
shifts on the nearly zero complex eigenvalues are used effectively, The average savings 
on test cases in this category is over 5 percent. 

In all the test cases presented in this paper, the stable eigenvalues from the two 
algorithms are exact to approximately 12 significant figures. Neither algorithm has 
shown consistently more accuracy than the other. Also, the unstable eigenvalues have 
been determined as accurately as possible by both algorithms. To determine the accu- 
racy of the unstable eigenvalues, one must check the accuracy of and the diag- 
onal elements of the resulting triangular A and B matrices, respectively. The 
and /3j are accurate up to a perturbation of order ||A||eQ and j|B| [6Q, respectively . 
For example, if the norms of the original matrices were of order unity and and j3j 
were of order 10 " then the eigenvalue Aj would have approximately one accurate digit 
even though Aj would be of order unity itself. To illustrate the point even further, test 
cases 1-4 and 1-5 for the 50 by 50 matrices are examples of eigenproblems which theoret- 
ically have all real eigenvalues, but both algorithms found some complex eigenvalues 
among the unstable eigenvalues. But, as expected, the imaginary parts of these eigen- 
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values had zero accuracy . Complex eigenvalues were the result of the nonpreserving- 
symmetry property which both algorithms possess. 

By looking at all the test cases, one can identify some definite trends concerning 
individual problems. They are as follows: 

(1) The more real eigenvalues a problem has, the more one can expect to save by 
using the combination shift QZ algorithm. Figure 1 gives a graphic view of this tendency. 
It shows a graph of the percentage of the QZ iteration time used by the combination shift 
QZ iteration plotted against the percentage of real eigenvalues . One standard deviation 
band about the average is also shown. 

(2) One can expect a larger savings on smaller matrices by using the combination 
shift QZ than on the larger matrices. This tendency is noted in practically all the test 
cases presented. Figure 2 gives a graphic view of the average and one standard devia- 
tion band of this trend. Also, as noted earlier, the standard deviation is larger for small 
matrices because of the greater effect of one iteration and the ” nonsettling" of the eigen- 



Figure 1.- Iteration time comparison by percentage of real eigenvalues. 
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values; that is, because of the lack of enough previous iterations, the eigenvalues have not 
settled down to a specified location or order. 

(3) The standard deviation tends to be larger for the cases which have a more equal 
distribution of real and complex eigenvalues. This difference may be attributed to a 
higher rate of reordering real and complex eigenvalues so that both shifts in the 

QZ algorithm are used effectively. 

(4) If the real eigenvalues have smaller magnitudes than the complex eigenvalues, 
one can expect a greater savings by using the combination shift QZ since the smaller 
eigenvalues tend to be found first. This tendency supports the comments presented in 
the previous section on the theoretical comparison. 

The results presented in the tables correspond to eigenvalue computation only . In 
the test cases which have been tried, approximately the same percentages resulted when 
the eigenvector calculation was added. The eigenvectors were calculated as they were in 
Moler and Stewart’s algorithm and not by the preferred method suggested by Kaufman. 
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As does the QZ algorithm, the combination shift QZ algorithm handles the "ill- 
disposed" problem when det(A - XB) vanishes identically, that is, when any X can be 
considered as an eigenvalue. This case appears with an "essentially" zero diagonal ele- 
ment on the final triangular matrices at the same relative location. 

CONCLUDING REMARKS 

The algorithm presented in this paper, called the combination shift QZ algorithm, 
solves the generalized eigenvalue problem. It should be used when both matrices are 
singular or nearly singular and tests indicate it is particularly effective on eigenprob- 
lems which have a large percentage of real eigenvalues. Based on the results presented 
in this paper, it should be preferred over existing algorithms which attempt to solve this 
class of eigenproblems. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 7, 1973. 
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APPENDIX 


TEST CASES 


Test cases were generated so that they fell into six basic categories. These cate- 
gories were pairs of matrices which in the generalized eigenproblem sense have: 

Category I - 100 percent real eigenvalues 

Category II - 80 percent real eigenvalues 

Category III - 60 percent real eigenvalues 

Category IV - 40 percent real eigenvalues 

Category V — 20 percent real eigenvalues 

Category VI - 0 percent real eigenvalues 

Within each category, six matrices were generated as a function of the matrix size N, 

To help present the details on the generation of the test cases, several matrices 
need to be defined. First, three orthogonal matrices and a tridiagonal matrix are as 
follows: 


i] 7J- 

U symmetric matrix whose ij element is sin - 

T matrix with diagonal equal to 10 and both Superdiagonals and subdiagonals 

equal to 4 


V orthogonal eigenvector matrix of T 


P orthogonal eigenvector matrix of a symmetric random number matrix with 

random numbers uniformly distributed in the interval |^-5., 5.j 

Next, several diagonal and block diagonal matrices, which are functions of the variables 
indicated and which define the eigenvalues, are defined as follows: 

Dj(N) = diag{l, 3, 5, . . ., N - 1; 10“ lO-H, 10-12; 

-N + 6, -N + 8, -N + 10, . . ., -4, -2} 

DgCN) = diag{l, 2, 3, . . ., N - 2; lO-H, 10-12} 
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APPENDIX - Gontinued 


Dg(N) = diagjl, 2, 3; + l)l0-12, (E - l)l0-12, . . (4)10-12; 


1, 2, 3, . . vf- 1 


D4(N) = diagjs, 5, 7; 4, 5, 6, . . |+ 1; 10-12, io-12, io-12, . . 10-12^ 


r 


Dg(N,k) = diag/ 


0 

-N 


N 

0 


0 N - 2 

-(N - 2) 0 


0 N - 4 

-(N - 4) 0 


0 
kN 


M ^2 

5 


+ 2 


. ^ M 1 Mi 9 
’ 5 ’ 5 ■ ’ 5 ■ ’ ■ 


3; 10-10, lo-liy 


Dg(N) =diag{l0-ll, 10-11; 1, 1, 1, . . ., l} 


r 


D;^(N,k) = diag/N, N - 1, N - 2 


M + i - 

■> C ^ 


kN 

5 


5 


0 
kN 


f -21 

5 


Mi_4 





0 

4 


, • . . 



5 


-4 

0 







0 

• 10-11 


10-11 

0 


D8(N) =diag{l0-ll, 10-10; 1, 1, 1, . . l} 
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Dg(N,k,Z,m,q) = diag/lO-H, 10-10; 3 ^ 4, 5, . . k; 


0 k+ 1 

-(k + 1) 0 


0 

-(k + 3) 


V 


0 

k + 5~ 


0 

Z - 1 



9 • • 



-(k + 5) 

0 


-{1 - 1) 

0 


; Z + 1, Z + 2, Z + 3, 


0 

m + 1 


0 

m + 3 


0 

m + 5 

-(m + 1) 

0 

9 

-(m + 3) 

0 


-(m + 5) 

0 


__ 



0 

q - 1 

; q + 1, q + 2, q + 3, . . N> 

-(q - 1) 

0 






Dio(N) = diag{l, 1, 1, . . .,1; IQ-IO, lO-H) 


Djj(N,k,Z,m,q) = diag< 


0 

10-11 


0 

3 


0 

1 

5 


1 

0 

7 

0 

t— t 
I 

0 



0 


-3 

1 

0 

9 

-5 

I 

0 

9 

-7 

0 


0 k - 1 

(k - 1) 0 


; k + 1, k + 2, k + 3, • . Z; 


k + 3 
0 


(Equations continued on next page) 
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-U + 1) 


U + 3) 


il + 5) 


m - 1 


(m - 1) 


; m + 1, m + 2, m + 3, , . q; 


q + 3 


q + 1 


-(q+1) 0 ~(q + 3) 0 -(q + 5) 


q + 5 


N - 1 


-(N - 1) 


0 N 


-N 0 -(N - 2) 


N - 2 


N - 4 


(N - 4) 


0 4 


-4 0 -10 


10-11 


10-10 


1.000000001 


0 - 1.000000001 


- 1.0001 


1.0001 


(Equations continued on next page) 
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APPENDIX - Continued 


2. 

7. 


2. 

9. 







> 

-7. 

0 


-9. 

0. 



2 . 11 . 
- 11 . 0 




2. N - 1 

-(N - 1) 0 


n- 


D^ 4 (N,k) = diag<( 


0 

■(k + 5) 


k -1- 5 
0 


0 k-H 7 

-(k-i-7) 0 


0 

-(N - 1) 


0 

10-10 

> 

0 

10-10 


0 

10-10 

-10-10 

0 


-10-10 

0 


-10-10 

0 









0 

-10-10 

10-10^ 

i 


0 

k-f 1 

? 

0 

k-f 3 

0 


-(k+1) 

0 

1 

-(k-f 3) 

0 


1 r 



-1 

1 

r 





N - 1 
0 


j 


Dl 5 (N,k) = diag( l, 1, 1^. ■ 1 ; IQ-lO, 10-10 , iq-10, . . iq-10) 

k -I- 2 components 


Some of the block diagonal matrices are constructed from block diagonal submatrices of 
dimension 10. They are as follows: 


r 


D^(N) = diag/lO-ll, 10-10, N - 2, 


0 

■(N - 3) 


N - 3 
0 




; N - 5, N - 6, N - 7, N - 8, N - 9\ 
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APPENDIX - Continued 


D 2 (N,k) = diag/N - 10k, N - 10k - 1, N - 10k - 2, 


N - 10k - 3 


-(N - 10k - 3) 


N - 10k - 5, N - 10k - 6, N - 10k - 7, N - 10k - 8, N - 10k - 9 


0 7 

Dg = diag/10, 9, 8, , 5, 4, 3, lO’lO, IQ-H 

-7 0 


0 N - 2 

D^CN) = diag/10-11, lQ-10, , N - 4, N ^ 5, 

1 -(N - 2) 0 


N - 6 


-(N - 6) 


, N - 8, N - 9 


Dg(N,k) = diag/N - 10k, N - 10k - 1, 


N - 10k ^ 2 


(N - 10k - 2) 


N ^ 10k - 4, N - 10k - 5, 


N - 10k - 6 


(N - 10k - 6) 


N - 10k - 8, N - 10k - 9 
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APPENDIX - Continued 




10-11 

N 


0 

N - 4 

D^(N) = diag< 

< 

-N 

10-10 

, N - 2, N - 3, 

-(N - 4) 

0 


■>> 

0 N - 8 

> 

-(N - 8) 0 

J 
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APPENDIX — Continued 


D^q(N) = diag< 


(r- 



p- 

- 



-] 

10-11 

N 

, N - 2, 

0 

N - 3 

> 

0 

N - 5 

-N 

10-10 


-(N - 3) 

0 


-(N - 5) 

0 


-J 


L, 






N - 7, 


0 

•(N - 8) 


N - 8 
0 




Dj^(N,k) = diag^ 



-| 


— 


0 

N - 10k 

, N - 10k - 2 , 

0 

N - 10k - 3 

-(N - 10k) 

0 


-(N - 10k - 3) 

0 


- 


- 

- 





_ 


0 

N - 10 k - 5 

, N - 10k - 7, 

0 

N - 10k - 8 

> 

~(N - 10k - 5) 

0 


- (N - 10k - 8) 

0 




- 



c- 


Di 2 = diag<[ 


0 

-10 


10 

0 


, 8 , 


0 

-7 


0 

-5 


, 3, 




10-10 

-2 


2 

10-11 


The test cases can now be defined as 


I-l: A = UT D^(N) U 

B = uT DgCN) U 

1-2: A = pTd^(N)P 

B = PT D2(N) P 

1-3: A = vTDj(N)V 

B = VT DgCN) V 


N = 10, 20, 30, 40, 50 
N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 
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APPENDIX - Continued 


1-4: 

1-5; 

1 - 6 ; 


II- 1: 


II- 2: 


n-3: 


A = UT DgCN) u) 
B = UT D^CN) u) 

A = DgCN) V) 
B = yT D4(N) V) 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


A = Symmetric random number matrix with random numbers uniformly 
distributed in the interval [^10., lOj 


B = CTc where the first N-3 rows of C are random numbers 
uniformly distributed in the interval [;5., 5^ and the last three rows 
are linear combinations of the preceding rows with a perturbation on 
the order of 10" added to each element in these rows. 


A = T D5(N,4) U 
B = T Dg(N) U 

A = T D^(N,1) U 
B = T Dg(N) U 

N = 10 


N = 20 


N = 30 


N = 40 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


(a = T DjClO) U 
(B = T Dio(lO) U 


A = T diag|D4(20), 03(20, 1)J U 
B - T Dio(20) U 


(a = T diag|D^(30) , 03(30, 1), 03(30 , 2)| U 
|b = T D4q(30) U 


^A = T diag|Di(40), 03(40, 1), 03(40, 2), 


0 , 


j(40, 3)} 


U 


B = T 0^0(40) U 


N = 50 


! A = T diag|Dj^(50), 03(50, 1), 03(50, 2), 
03(50, 3), 03(50, 4)| U 
B = T Dio(50) U 
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APPENDIX - Continued 


II-4: 


n-5: 


n-6: 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


Let 


(a = T Dg U 
(B = T Dg(lO) U 

(a = T diag|D^2(20, 0), Dg| U 
I B = T Dg(20) U 

(a = T diag{D2(30, 0), 03(30, 1), Dg| U 
)b = T Dg(30) U 

(a = T diag|D2(40, 0), 03(40, 1), ^3(40, 2), Dg| U 
Ib = T 03(40) U 

! A = T diag{02(5O, 0), ^3(50, 1), ^3(50, 2), 

02(50, 3), Dg| U 
B = T Dg(50) U 

(a = T Dg(10, 4, 6 , 10, -) U 
(b = T Diq(10) U 

(a = T Dg(20, 5, 9, 20, -) U 
(b = T D^q(20) U 


(A = T 09(30, 6 , 12, 30, -) U 
|b = T Oio(30) U 

(a = T Dg(40, 4, 10, 18, 20) U 
(B = T D jq(40) U 

(a = T 09(50 , 5, 11, 21, 23) U 
(B = T Diq(50) U 

a’ = A of test case II- 5 

B' = B of test case II- 5 

A = U T-1 A' U-1 t) 

) N = 10, 20, 30, 40, 50 

B = U T-1 B’ U-1 T) 
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APPENDIX - Continued 


in-1: 

III- 2: 
III- 3: 


m-4; 


A = T D5(N,3) U 
B = T Dg(N) U 

A = T D7(N,2) U 
B = T Dg(N) U 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


N = 10 


N = 20 


N s 30 


N = 10, 20, 30, 40, 50 

N = 10, 20, 30, 40, 50 

(a = T D^(IO) U 
(B = T Diq(10) U 

(a = T diag|D4(20), 05(20, 1)| U 
]b = T D4o(20) U 

(a = T diag|D4(30), 05(30, 1), 05(30 , 2)| U 
|b = T D^qOO) U 

i A = T diag|D4(40), D-5(40, 1), 05(40,2), 
1X5(40, 3)| U 

B = T Diq(40) U 

i A = T diag{D4(50), D'5(50, 1), 05(50, 2), 
05(50, 3), 05(50, 4)| U 

B = T 0^0(50) U 

( a = T Og U 
(b = T Dg(lO) U 

(a = T diag|o5(20, 0), Dg} U 
Ib = T Dg(20) U 

(a = T diag|D5(30, 0), Og(30, 1), Ogj U 
|b = T Dg(30) U 
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APPENDIX - Continued 


m- 5 : 


m-6: 


IV-1: 

IV-2; 

IV-3: 


N = 40 


N = 50 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


Let 


A = T diag|D5(40, 0), 05(40, 1), 05(40, 2), Dg| U 
B = T Dg(40) U 


A = T diag{o5(50, 0), 05 ( 50 , 1), 05 ( 50 , 2), 

05(50, 3 ), C-g j U 


B = T Dg(50) U 


|A = T Dg(10, 5, 9, 10, -) U 
[b = T Dio(lO) U 


JA = T 09(20, 7, 15, 20, -) U 

(B = T Dio(20) U 

(a = T Dg(30, 5, 13, 21, 25) U 

(B = T Dio(30) U 

(a = T 09(40, 8, 16, 24, 32) U 

(B = T Dio(40) U 


|A = T 09(50 , 7, 15, 27, 39) U 
[b = T Oio(50) U 


A' = A of test case 111-5 
B' = B of test case IH-S 
A = U T-1 A' U-1 t) 

B = U T-1 B’ U-1 T) 


N = 10, 20, 30, 40, 50 


N = 10 


A = T D5(N,2) U 
B = T Dg(N) U 

A = T D7(N,3) U 
B = T Dg(N) U 

A = T O'7(10) U 
B = T Diq(10) U 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 
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APPENDIX - Continued 


IV-4: 


IV- 5: 


N = 20 


N = 30 


N = 40 


N = 50 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


N = 10 


= T diag{D7(20), DgCBO, 1)} U 
B = T D U 

'a = T diag|D7(30), UgCSO, 1), ITgCSO, 2)| U 
B = TDjq(30)U 

I A ^ T diag{D7(40), D8(40, 1), Dg(40, 2), S'8(40, 3)| U 
I^B = T Dio(40) U 

A = T diag|D7(50), D8(50, 1), DgCSO, 2), D'8(50, 3), 
D8(50, 4)| U 
B = T Dio(50) U 

jA - T Dg U 
(B = T DgClO) U 


A = Tdiag|D8(20,0),D9|U 
'b s T DgC^O) U 

I A T diag|58(30, 0), UgCSO, 1), Dg| U 
B - T DgCSO) U 

A = T diag|D8(40, 0), Dg(40, 1), D8(40, 2), Dg^ U 
'B = TDg(40)U 

V 

/ 

A = T diag|D8(50, 0), D8(50, 1), Dg(50, 2), 

D8(50, 3), Dg|u 
B = T DgCSO) U 

(A = T DgClO, 3, 9, 10, -) U 
|B = T U 
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APPENDIX — Continued 


IV- 6: 

V- 1: 
V-2: 
V-3: 


(A = T Dq(20, 3, 9, 13, 19) U 
N = 20 \ 

(B = T Dio(20) U 

(a = T Dq(30, 6, 16, 20, 28) U 
N = 30 { ^ 

(b = T Dio(30) U 


N = 40 


N = 50 


(A = T Dg(40, 6, 18, 22, 34) U 
(B = T Dio(40) U 
(A = T Dg(50, 10, 24, 30, 46) U 
(b = T D^gCSO) U 


Let a' = A of test case IV- 5 
b' = B of test case IV- 5 
A - U T -1 A’ U-1 T 
B = U T-1 B’ U-1 T 


N = 10, 20, 30, 40, 50 


A ~ T D5(N,1) U 
B = T Dg(N) U 

A = T D,y(N,4) U 
B = T Dg(N) U 

(A=TDio(10)U 
N = 10 { 

B = T Diq(10) U 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


N = 20 


N = 30 


N = 40 


(a = T diag|D^Q(20), Dii(20, 1)^ U 
|b = T Dio(20) U 

j A = T diag|DjQ(30), Dn(30, 1), Dn(30, 2)| U 
|b = T Diq(30) U 

(a = T diag|Dio(40), 1Jii(40, 1), nii(40, 2), ITii(40, 3)| U 
(B = TDio(40) U 
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APPENDIX - Continued 


V-4: 


V-5: 


N = 50 


N = 10 


N = 20 


N = 30 


N = 40 


N = 50 


N = 10 


N = 20 


N = 30 


N = 40 


(a = T diag|Dio(50), 1), Di^(50, 2), E 

< Dii(50, 4)| U 

(^B = T Dio(50) U 

|A = T Di 2 U 
(b = T Dg(lO) U 

f A = T diag|Dn(20, 0), U 
|b = T Dg(20) U 

(a = T diag|Dij^^(30, 0), Pii(30, 1), D^2} U 
| b = T Dg(30) U 

|a = T diag|Dii(40, 0), D-ii(40, 1), 2) 

/ B = T Dg(40) U 

i A = T diag|Dii(50, 0), %(50, 1), 2) 

Dii(50, 3), D12} U 
B = T Dg(50) U 


(A = T D^^(10, 6, 8, 10, -) U 
|b = tDio(io)u 
(A = T Dij(20, 6, 10, 20, -) U 
( B = T D^o(20) U 

I A = T Dii(30, 8, 14, 30, -) U 
[b = T D^q(30) U 

|A = T Dii(40, 6, 10, 20, 24) U 
(B = T Dio(40) U 


jj^(50, 3), 


.*^12}U 
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APPENDIX - Concluded 


N = 50 


V-6: 


VI- 1: 


VI- 2: 


VI-3: 


VI-4: 


VI- 5: 


VI- 6: 


Let 


Let 


Let 


(A = T D^j(50, 10, 16, 24, 28) U 
[b = T Dio(50) U 

A* = A of test case V-5 
B’ = B of test case V-5 
A = U T-1 A’ U-1 T 
B = U T-1 B’ U-1 T 

A = T Di 2(N) U j 
B = T D0(N) U 

A = U Di2(N) T | 

B = U Dg(N) T 

A = T Di3(N) U | 

B = T Dio(N) u) 

A = UDi3(N) 

B = U Diq(N) t| 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


N = 10, 20, 30, 40, 50 


m = 2 


where jr] means the greatest integer not 
greater than r 
A = T Dj4(N,m) U 
B = T Di5(N,m) U 
m = m of test case VI- 5 
A = U Dj4(N,m) T 
B = U Di5(N,m) T 
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TABLE I.- ITERATION TIME COMPARISON FOR CASES WITH 
ALL REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 

Row 

standard 


N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

deviation 

I-l 

70.5 

61.9 

58.8 

68.0 

65.2 

64.9 

3.9 

1-2 

44.9 

52.1 

67.5 

65.3 

65.6 

59.1 

8.8 

1-3 

55.8 

55.4 

66.2 

68.0 

65.0 

62.1 

5.1 

1-4 

62.8 

64.7 

62.8 

67.2 

68.3 

65.2 

2.3 

1-5 

62.8 

66.1 

61.9 

66.5 

70.1 

65.5 

2.4 

1-6 

69.0 

72.1 

73.5 

69.8 

73.8 

71.6 i 

3.1 

Column average . . . 
Column standard 

61.0 

62.1 

65.1 

67.5 

68.0 



deviation 

8.4 

6.2 

4.9 

1.6 

3.2 




Average percentage for these cases 64.7 

Standard deviation for these cases 6.4 


Average percentage for these cases 64.7 

Standard deviation for these cases 6.4 




TABLE IL- ITERATION TIME COMPARISON FOR CASES WITH 
80 PERCENT REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 

Row 

standard 


N = 10 

N = 20 

N = 30 

N = 40 

N - 50 

deviation 

IM 

58.0 

70.7 

78.6 

77.6 

77.9 

72.6 

7.4 

II-2 

104.7 

97.7 

80.9 

88.1 

74.6 

89.2 

10.9 

II- 3 

60.5 

63.6 

69.8 

70.3 

76.5 

68.1 . 

6.1 

n-4 

60.5 

69.3 

71.9 

73 ;8 

76.8 

70.5 

5.0 

II- 5 

70.7 

62.0 

70.0 

78.0 

69.0 

69.9 

5.6 

II- 6 

65.7 

85.9 

85.6 

74,2 

85.2 

79.3 

8.3 

Column average . . . 
Column standard 

70.0 ' 

74.9 

I 

76.1 

77.0 

76.7 

] 


deviation ...... 

12,7 

12.6 

6.4 

5.6 

4.2 



Average percentage for these cases . , 

.... 

. . . . . 

> .... 


, . . 74.9 

Standard deviation for these cases . . . 

.... 

.... 

. . . . . 


. . , 10.7 
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TABLE III.- ITERATION TIME COMPARISON FOR CASES WITH 
60 PERCENT REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of — 

Row 

average 

Row 

standard 

deviation 

N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

m -1 

60.0 

71.2 

71.3 

82.0 

84.9 

73.9 

8.7 

m -2 

104.4 

94.2 

92.4 

77.7 

77.8 

89.3 

10.3 

m-3 

58.3 

65.9 

72.4 

66.5 

71.6 

66.9 

5.6 

m-4 

58.3 

65.9 

66.8 

61.8 

71.1 

64.8 

4.1 

m-5 

46.2 

70.5 

70.3 

73.3 

78.9 

67.8 

11.5 

m-6 

63.2 

49.5 

52.1 i 

65.0 

77.9 

61.5 

1 

10.4 

Column average . . . | 

65.1 

69.5 

70.9 1 

71.1 

77.0 



Column standard 








deviation 

18.2 

13.3 

11.7 

6.7 

5.2 




Average percentage for these cases . 70,7 

Standard deviation for these cases 12.7 


Average percentage for these cases . 70,7 

Standard deviation for these cases 12.7 
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TABLE IV.- ITERATION TIME COMPARISON FOR CASES WITH 
40 PERCENT REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 

— 

Row 

standard 

deviation 

N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

IV-1 

60.0 

60.1 

70.4 

84.4 

81.7 

71.3 

10.5 

IV- 2 

112.5 

118.8 

94.8 

94.7 

90.4 

102.2 

11.6 

IV-3 

66.7 

59.8 

73.2 

73.1 

80.4 

70.6 

7.3 

IV-4 

61.3 

73.5 

78.8 

75.8 

79.9 

73.9 

6.2 

IV-5 

1 

78.0 

84.7 

76.9 

82.8 

77.0 

79.9 

2.7 

IV-6 

101.6 

97.0 

89.0 

89.9 ; 

90.7 

93.6 

5.6 

Column average . . . 

80.0 

82.3 

80.5 

83.5 

83.4 



Column standard 








deviation 

20.3 

21.0 

8.8 

6.9 

4.4 




Average percentage for these eases 81.9 

Standard deviation for these cases 14.5 


Average percentage for these eases 81.9 

Standard deviation for these cases 14.5 
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TABLE V.- ITERATION TIME COMPARISON FOR CASES WITH 
20 PERCENT REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 

Row 

standard 


N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

deviation 

V-1 

81.4 

91.4 

79.2 

83.7 

87.4 

84.6 

4.7 

V-2 

100.0 

107.9 

104.2 

100.5 

100.4 

102.6 

3.1 

V-3 

96.8 

85.3 

82.4 

86.1 

84.5 

87.0 

5.4 

V-4 

91.8 

87.4 

86.0 

87.0 

81.5 

86.7 

4.2 

V-5 

88.9 

66.4 

84.3 

83.6 

94.9 

83.6 

9.7 

V-6 

41.7 

64.8 

83.8 

97.0 

98.3 

77.1 

21.5 

Column average ... 
Column standard 

83.4 

83.9 

86.7 

89.6 

91.2 



deviation ^ 

19.7 

14.6 

7.6 

7.3 

i 

6.7 




Average percentage for these cases 87.0 

Standard deviation for these cases 12.4 


Average percentage for these cases 87.0 

Standard deviation for these cases 12.4 
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TABLE VI.- ITERATION TIME COMPARISON FOR CASES WITH 
ALL COMPLEX EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 

Row 

standard 


N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

deviation 

VI- 1 

102.4 

101.1 

101.1 

100.1 

99.5 



VI-2 

82.9 

88.3 

92.5 

95.0 

88.3 

89.4 

4.1 

VI- 3 

85.0 

100.8 

94.8 

95.3 

96.2 

94.4 

5.5 

VI-4 

102.4 

89.5 

99.1 

95.4 

98.8 

97.0 

5.2 

VI- 5 

100.0 

101.1 

101.0 

100.9 

100.6 



VI- 6 

75.6 

91.8 

75.2 

87.1 

95.7 

85.1 

8.2 

Column average ... 
Column standard 

91.4 i 

95.4 

94.0 

95.6 

96.5 



deviation ...... 

10.5 

6.2 

8.4 

5.2 

4.4 




Average percentage for these cases 94.6 

Standard deviation for these cases 7.2 


Average percentage for these cases 94.6 

Standard deviation for these cases 7.2 
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TABLE VII.- OPERATION COUNT COMPARISON FOR CASES WITH 

ALL REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 


N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

I-l 

63.0 

58.9 

56.7 

67.0 

64.7 

62.1 

1-2 

39.6 

48.9 

65.3 

64.3 

65.3 

56.7 

1-3 

49.5 

52.2 

64.3 

67.0 

64.5 

59.5 

1-4 

56.1 

60.6 

59.8 

64.1 

66.4 

61.4 

1-5 

! 

56.1 

62.1 

58.7 

64.4 

68.2 

61.9 

1-6 

60.2 

68.5 

71.5 

68.4 

73.2 

68.4 

Column average .... 

54.1 ! 

58.5 

62.7 

65.9 

67.1 


Average percentage for these cases . . . 

• • ' • 



. . 61.7 




TABLE Vni.- ALGORITHM TIME COMPARISON FOR CASES WITH 

ALL REAL EIGENVALUES 


Test case 

Percent of QZ used by combination 
shift QZ for matrix size of - 

Row 

average 


N = 10 

N = 20 

N = 30 

N = 40 

N = 50 

I-l 

82.7 

77.6 

76.3 

82.5 

81.2 

80.1 

1-2 

65.0 

71.4 

81.8 

80.7 

81.5 

76.1 

1-3 

70.7 

73.9 

81.3 

82.5 

81.2 

77.9 

1-4 

79.5 

83.1 

84.7 

88.6 

89.5 

85,1 

1-5 

77.0 

84.0 

83.9 

88.2 

90.5 

84.7 

1-6 

79.5 

83.8 

84.9 

83.2 

85.3 

83.3 

Column average .... 

75.7 

78.9 

82.2 

84.3 

84.8 


Average percentage for these cases . . . 

.... 



. . 81.2 
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